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ABSTRACT 


Ship commanders and pilots make life or death decisions based on the set of 
information they have at their disposal at the instant a decision is made. One component 
of that information is the question of whether a radar contact is an enemy or a friend. 
Various systems exist which try to answer that question based on the characteristics of 
signals emitted or scattered from the contact. Any hesitation made by the operator as a 
result of not trusting the identification being presented by his system can be fatal. The 
goal is to maximize the quality of the information he receives in order to build the trust 
that when the system tells the operator the contact is an incoming friendly, he knows that 
it is. 

This thesis examines the technique of using the bispectrum of backscattered radar 
energy to identify a contact. Bispectra allow the examination of multiple scattering 
contributions to the backscattered energy. This technique is compared to one using radar 
range profiles. A library of sample radar signatures is built using computational radar 
cross section estimation tools and 3-D model aircraft. This library is the basis of a series 
of simulations with aircraft at multiple aspects and configurations to determine whether 
using the bispectrum enhances the performance of identification systems using range 
profiles. It is determined that a bispectrum method meets or exceeds the identification 


accuracy of the range profile method especially with high-bandwidth systems. 
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I. INTRODUCTION 


A. MOTIVATION 


The importance of having accurate information on the battlefield cannot be 
underestimated. Critical decisions must be made whether all the pertinent information is 
available or not. The decision to shoot down an incoming aircraft sometimes must be 
made without knowing exactly what that aircraft is. Inaccurate or unavailable 
identification information can therefore lead to instances where friendly or neutral forces 
are fired upon inadvertently. Improving the quality of information provided by the 
combat systems will enhance our warfighting capabilities allowing our warfighters to 


have a more accurate view of the battle space. 


1. Detect-to-Engage Sequence 


Every naval tactician is taught in their basic naval science class the detect-to- 
engage sequence. It is the sequence of steps through which every successful encounter 
with the enemy progresses, whether the sailor is on guard duty with a rifle or a Tactical 
Action Officer (TAO) on an Aegis Cruiser. I will look at this process from the point of 
view of the TAO. The first step in the process is detection. There are many ways a 
modern warship can detect a target, but a common method is via radar. Once the radar 
detects a target, it will enter it into its database and start to track it. Combining 
information from all available sensors, identification of the contact can be accomplished. 
Based on the current rules of engagement, a threat assessment of the contact is 
performed. A contact deemed hostile will be paired with a particular weapon or platform 
to commence the engagement. Following the engagement, a damage assessment must be 


performed to determine if an additional engagement is necessary. 


Each of these steps is not always a formalized process and the particular order can 
be rearranged to a certain extent. Our advancing weapons technology has made the 
engagement an increasingly successful part of the process, and modern radars provide a 
well over the horizon look at the airspace. The most troubling part lies with the 
identification and threat assessment of the contact. There is still an element of 


1 


subjectivity with the threat assessment; it can be based on experience and intuition and is 
very situation dependent. Many different techniques exist to identify contacts, but none 
are without flaw. Also, an enemy’s efforts in deception can change characteristics that 
would otherwise reveal their identity. These elements of the sequence are big parts of 
what keeps the process from being fully automated. Improving the reliability of the 


identification phase will also enhance the assessment of the threat. 


2. The Price of Failure 


There have been a number of incidents in modern American history which 
illustrate the consequences of having inaccurate identification information. These cases 
not only lead to the destruction of millions of dollars in property, but also lead to the 


deaths of many innocent civilians or friendly military forces. 


1988 USS Vincennes incident — On June 2, 1988, the Aegis cruiser USS 
Vincennes shot down Iran Air Flight 655. The Vincennes and the frigate USS 
Elmer Montgomery were engaged in a surface gun battle with Iranian gunboats. 
The Iranian Airbus A300 took off from Bandar Abbas and was misidentified by 
the Vincennes as an F14 fighter aircraft. Despite flying in a commercial flight 
path, the aircraft was labeled hostile and was shot down by two standard missiles 


from the Vincennes. All 290 persons aboard were killed [1]. 


1994 Iraqi no-fly zone — On April 14, 1994, two Air Force F-15 fighters shot 
down two American MH-60 Blackhawk helicopters in the no-fly zone in northern 
Iraq. The helicopters were not identifying themselves properly via their IFF 
system, (see Chapter IV) and they were misidentified electronically and visually 
by the fighters and surveillance aircraft. 26 military and civilian personnel lost 


their lives [2]. 


2003 Downing of British Tornado — On March 23, 2003, a U.S. Patriot missile 
battery shot down a British GR4 Tornado aircraft after misidentifying it as an 
incoming missile. The GR4 was returning from a mission in Iraq. The Patriot 
missiles can operate in mostly autonomous modes where automatic identifications 


play a major role in the engagement [3]. Both pilot and navigator were killed. 
2 


B. ORGANIZATION 
1, Chapter Descriptions 


This thesis examines a way to improve the accuracy of radar contact 
identification. Chapter II begins with an overview of basic radar concepts. Included is a 
brief discussion of different signal forms and high resolution techniques which enable the 
later technologies. Chapter III deals with how the radar waves actually interact with the 
aircraft and scatter back to the receiver. The different contributions to the radar 
backscatter are examined with a simple example. Chapter IV gives an overview of 
identification methods of radar contacts. Both cooperative and non-cooperative 
classification methods are examined including the shortcomings of the current methods. 
In addition the technique of using the bispectrum of the radar signal to improve the 
accuracy of classifications is described. A method of how to calculate the bispectrum 
and an example using simple shapes is shown. Chapter V describes the simulation that 
was performed to test the hypothesis that the use of the bispectrum will improve 
classification accuracy. A discussion of the creation of a library of aircraft signatures and 
a selection of sample test cases is included. Chapter VI provides results from the 
simulation. Additional observations using the library of signatures are made. 
Conclusions follow in Chapter VII to summarize the findings and to recommend further 


examination. 


2 Thesis Objective 


A library of radar backscattered data from various models of military and civilian 
aircraft is desired to test various classification methods. This library will be accessible 
for future work beyond this thesis. A comparison of radar range profile based 
classification methods versus a method based on the bispectrum of the responses is 
performed. The goal is to improve the reliability of automated radar based classification 


techniques by demonstrating a promising method to enhance current methods. 
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I. OVERVIEW OF RADAR PRINCIPLES 


A basic understanding of how a radar system works is helpful in understanding 
the goals of this thesis. The word radar began as an acronym for Radio Detection and 
Ranging but has become so pervasive that it is now a word by itself. Radar is a system 
which determines a contact’s range, bearing, velocity, and other information by emitting 
electromagnetic energy and measuring the energy that reflects from the contact. Radar 
signals can be based on pulses which alternate between transmitting and receiving 
energy, or they can be continuous wave (CW) systems which constantly transmit and 


receive simultaneously. 


The wide scale use of radar began in World War II. A well-known system was 
the British Chain Home radar which allowed advance notice of attacking German fighters 
and bombers [4]. British fighters were able to be vectored directly to where the incoming 
planes were, therefore greatly enhancing their effectiveness. The United States was also 
working on its own radar at the time. The invention of the magnetron to deliver high 
power, coherent microwaves, and the duplexer switch enabled construction of a pulsed 
radar which was much smaller than the original lower frequency radars [5]. Higher 
frequency radars allow the use of smaller antennas which make the system portable and 
more easily mounted on mobile platforms. Radar technology has drastically improved 
over the years with the advancement of signal processing techniques. Also, complicated 
antenna designs allow energy to be directed to desired locations without mechanically 


moving the antenna. 


A. BASIC HARDWARE SETUP 


Because of the spreading of electromagnetic waves, the reflected and received 
energy after the round trip is multiple orders of magnitude less than the transmitted 
energy. The receiver therefore must be a very sensitive piece of equipment. A duplexer 
allows it to be connected to the same antenna of the transmitter without being 
overwhelmed by the transmitted energy. A basic layout of a pulsed radar is depicted in 


Figure 1. Waveforms from the pulse generator allow varying levels of resolution and 
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will be discussed in the next section. Heterodyne mixing is done before transmitting and 
after receiving to step up an intermediate frequency to a higher frequency. This is done 
because a higher frequency transmitted may have better transmission and reflection 
characteristics, but the lower intermediate frequency is easier to manage for the electronic 
signal processing components. Various displays have been developed over the years 
beginning with phosphor scopes, and now high-tech multifunction computer displays 
found on modern warships. Filters and low noise amplifiers ensure the detector receives 


the right amount of power for detection without saturating or damaging components. 
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Figure 1. Simple Pulsed Radar [6] 


B. ANTENNA DESIGN 


After the energy is produced by the system, the antenna is what transfers that 
energy into the environment. The particular antenna design is motivated by the intended 
use of the system. These antennas could be simple dipoles to be used as a 360-degree 
proximity warning device for ships similar to Christian Huelsmeyer’s ship avoidance 
system in 1904 [5]. Or the antenna could be a complex, electronically steered, phased 
array antenna found on modern warships and aircraft generating a narrow beam. 
Traditional radar antennas consist of a feed horn at the end of a waveguide which directs 
energy to a parabolic dish. The shape of the parabolic dish largely determines the beam 
pattern formed. The electrical current distribution excited on the dish creates the 


radiation pattern and in general, the far-field beam pattern is the 2-dimensional Fourier 
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Transform of the aperture shape [7]. Because of the Fourier relationship between spatial 
frequency and length, wide dimensions provide narrow beams and narrow dimensions 
provide wider beams. Also, by adding a windowing function to the energy directed at the 
aperture, the sidelobes of the signal can be reduced, albeit at the cost of a wider main 
lobe. By grouping multiple elements together in an array, areas of constructive and 
destructive interference create a beam pattern that can be steered by changing the phase 
difference of the transmitted waveform between elements. A sample beam pattern 


showing sidelobes is displayed in Figure 2. 
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Figure 2. Sample beam pattern for rectangular aperture. [6] 


C. THE RADAR EQUATION 


The radar equation estimates the theoretical maximum range at which a radar can 


detect a target. If we take the power density of the transmitted antenna to be P. then, 


taking into account the gain of the antenna, G, and spherical spreading of the waves, the 
power density at the target is given by [5]: 
P.G 
ane = rey 


This is the power available to be scattered from the target. The amount of energy 


(1) 


scattered from the target is given in terms of its radar cross section, o , with dimensions 
of area. Assuming this reflected energy expands omnidirectionally, it will be further 
reduced by a 1/R’ factor. The antenna can only collect the fraction of power that hits the 
dish, so this effective area, Acs, is then multiplied to give the total power received by the 


radar: 


ae cas 2) 
Based on the electronics and signal-to-noise ratio (SNR) of the radar, if Pmin is the 
minimum detectable power, then rearranging (2) to solve for R, we get the maximum 
theoretical range for the radar to be: 

P,GAj-o ‘ 


R... =| ———_ 3 
max (47) P (3) 


This range is of course only theoretical. We have not taken into account any 
atmospheric effects like attenuation or the effect of the shape of the target on scattering. 
Also, care must be taken to avoid ambiguous range in the returns. If a new pulse is sent 
out before the round trip time of the maximum range, then the system may consider the 


return a result of the new pulse and provide a much closer range than exists. Techniques 


have been developed to avoid this problem [8]. 


D. RESOLUTION 


Good resolution allows a radar system to distinguish one large contact from many 
smaller contacts (i.e.,. planes flying in formation). If the resolution is high enough, then 
characteristics on individual planes can be distinguished. The two important resolutions 


in radar are bearing resolution and range resolution. 


1. Bearing Resolution 


Bearing resolution is a measure of the minimum discernable angular separation of 
contacts. It is determined by three factors. First, the beam width and therefore the 
aperture play a role. Obviously an omnidirectional antenna will provide no bearing 
resolution at all. The narrower the beam, the better the resolution. Second and third are 
the rotation rate of the antenna and pulse repetition frequency (PRF). These two must be 
optimized together. If the antenna rotates more than a beam width between pulses then 
the radar system will not have complete coverage of the area. It may completely miss a 
target between pulses. Greater desired ranges mean a slower PRF and therefore a slower 


rotation rate for complete coverage. 


2 Range Resolution 


What enables the techniques in this thesis is good range resolution. This is the 
minimum radial distance between 2 objects that can be differentiated by a single radar 
pulse. The range resolution will be half the width of the pulse if we think of a simple 


pulse depicted by the rectangular function in Figure 3. 
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Figure 3. Range resolution using simplified waveform 


Much more clever schemes have been derived based on pulse compression and 
correlation receivers. A common way to achieve greater range resolution is by using 
more complex waveforms [9]. For instance, instead of using a simple sign wave pulse, a 
pulse containing changing frequencies can be used. It is called an FM chirp and is 


depicted in Figure 4. 
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Figure 4. FM chirp representation. [6] 
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When this non-symmetrical waveform scatters from contacts, the received energy 
is then correlated with a time and frequency shifted copy of the transmitted waveform. 
Much closer objects are able to be resolved in this manner. The minimum resolvable 
separation, Or , using this scheme is given by 


c 
jer (4) 
It is clear that in order to increase the range resolution of these radars, more bandwidth is 
needed. It is common for bandwidth to be referred to in terms of fractional bandwidth; 
that is, the fraction of the total bandwidth over the carrier frequency. For instance, a 
system which operates around 10GHz with a 1GHz bandwidth will have a fractional 


bandwidth of 10%. 


This has been a very brief overview of radar fundamentals and many more 
considerations exist. The important information applicable to the understanding of the 
following material has been highlighted. I have discussed how energy is transmitted and 
received and next we will discuss just what happens when that energy encounters a 


target, and how that energy is reflected. 
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Hl. ELECTROMAGNETIC SCATTERING 


I’ve discussed transmitting energy and waveforms, but what is this energy and 
how does it propagate? In order to discuss how electromagnetic waves scatter from radar 
targets, we must first analyze how they get there in the first place. As in most endeavors 
in electromagnetism, the answer lies in Maxwell’s Equations. How do these equations 
yield the wave equation, and what are the possible solutions? Next, how do these waves 
interact with the structure of the targets and reflect waves back to the antenna? Lastly, 


what information can be inferred about the target from the received energy? 


A. ELECTROMAGNETIC WAVES 


The relationship between fields and charges is described by Maxwell’s Equations: 


V-D=p, V-E=0 
V-B=0 V-B=-0 
vxpu 2B Verso 

Ot Ot 
Sa yes SOD) OE 
VxH= A VxB Hoe, 
! 6t a) 


The second set of equations are reached as a result of making a simplifying 
assumption that we are dealing with a source-free, non-conductive, linear, homogenous, 
and isotropic medium [10]. These equations are the fundamental equations of 
electromagnetism and allow us to understand an electromagnetic wave as a time-varying 
disturbance of the electric and magnetic field which propagates at the speed of light 
through a medium in a direction perpendicular to each of the fields. A typical 


representation of the wave is depicted in Figure 5. 
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Figure 5. Representation of an Electromagnetic Wave. [11] 


If we take the curl of Maxwell’s third equation, apply a vector identity, and make 
substitutions we obtain the vector wave equation 
OE 

ot 


If we assume that solutions to the wave equation (and Maxwell’s equations for that 


VE-LE, =0 (5) 

matter) have harmonic time dependence, then the solution can be written in the form 
E(x,t) = E(xe (6) 

Plugging equation (6) into (5), the wave equation becomes the reduced wave equation, 


also known as the Helmholtz equation 





(Vv +k?) E(x) =0 where k = @ and c = ie " 


Further simplification of the solution for the wave (6) can be made if we assume 
that the target is very far away from the source. This is a very good assumption for radar 
since the target is normally well outside visual range. In this case, the radius of the 
spherically expanding wave is so great that from the targets perspective, the incoming 
wave is planar. This approximation allows us to say that the wave has no dependence on 
coordinates perpendicular to the direction of propagation. The solution can then be 
written with harmonic spatial and temporal frequencies. 
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EG@ijene (8) 
This form of a wave is a theoretical abstraction since a perfect plane wave does not exist 
in nature. For a local frame, far away from the source, it is a reasonable assumption and 
it will be used to model the incident wave for the scattering discussion in the next section 


[12]. 


B. SCATTERING 


An electromagnetic wave incident on a conducting surface is a complicated 
process. Even if we consider a perfectly electrically conducting surface (PEC), to 
simplify the boundary conditions, the analysis is still mathematically intensive. Although 
many times it is thought of as the wave ‘bouncing’ off, this is a misleading perception. In 
introductory optics, we draw rays hitting mirrors and being redirected as if they’re 
billiard balls on a pool table. Whereas this is a reasonable approximation to make in a 
high frequency regime, at radar frequencies it does not tell the full story. Upon closer 
examination, phenomenon such as evanescent waves, skin depth, induced current 
distributions, and interference are observed. When an electric field reaches a conductor, 
charges in the conductor are accelerated by the field. We consider a PEC whose charges 
experience no resistance in motion and can therefore respond immediately to the incident 
field. These moving charges are a current distribution and therefore also create radiation 
in the form of a scattered field. Together, the incident field and the scattered field must 
always satisfy Maxwell’s Equations and in turn, the wave equation. The time- 
independent solutions to the scalar wave equation in spherical coordinates will have the 
form 

T17' (r,0,¢) =P (kr) Y," (0,4) (9) 
where Z refers to a Bessel Function and Y to a spherical harmonic. Finding these 


solutions is very difficult for even the most trivial shapes [13]. 


Consider the plane wave of equation (8) incident on a PEC sphere. It can be 
shown that the scattered energy versus wavelength behaves as in Figure 6. In the 


Rayleigh region, the sphere is much smaller than the wavelength and the scattered energy 
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is proportional to @*. In the Mie region, the sphere and the wavelength are of the same 
order of magnitude giving rise to resonances as the electrons are accelerated producing 
constructive and destructive interference. In the optical region, the object is much larger 
than the wavelength and it becomes more valid to consider the incident wave as rays 


bouncing off the surface. 
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Figure 6. Scattering of a Plane Wave on a PEC Sphere [14]. 


Since developing analytical solutions to complex scattering geometries is beyond 
our capabilities, many techniques have been devised to enable automation via a computer 
to computationally approximate scattered fields from arbitrary objects. Techniques used 
in the generation of the datasets for this thesis are physical optics (PO) [15], the physical 
theory of diffraction (PTD) [16], and multiple interactions using shooting and bouncing 
rays (SBR) [17]. Together, these methods try to approximate scattered fields in a 


relatively high frequency regime. 


1. Physical Optics 


The physical optics method calculates the scattered field by relating the incident 
field to the induced surface currents. This method only considers the initial bounce of the 
incident wave and not contributions from energy multiply scattered from other parts of 
the object. This method also does not take into account fringe effects. Since all of the 
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aircraft models used in the simulations are assumed to be PEC objects, this discussion 
will also omit the effects of dielectric coatings. If we make what is called the weak 
scatterer approximation (also known as the first Born approximation), ignoring effects 
from multiple interactions, then we can write the current distribution as a result of the 


incident field as 


bl! 


(Xt) © J po (%,t) = 20x H inc (%, 1) (10) 
Using this current density, we can calculate the vector potential A and, in turn, the 
scattered magnetic field. After much algebra, the equation becomes 


*7E7 i(2kR-ot 
—ikH eC") 


2nR ang Re as! ny) 


A seat (Rtk) = 


R is the direction from the radar to the target and the integral is over the illuminated part 
of the object. In practice, the target’s smooth shape is discretized into hundreds of 
thousands of small triangles called facets as illustrated in Figure 7. The computer 
program calculates this integral over each of these flat surfaces and sums the 
contributions. Representing smooth objects with triangles, of course, introduces error 
into the calculation, but if modeled carefully with small triangles, this error can be 
minimized. A ray tracing algorithm is used to ensure only those illuminated facets are 


included in the calculation. 





(a) Wireframe (b) Flat Shaded 


Figure 7. Facetization of a Surface [18] 
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2. Physical Theory of Diffraction 


Physical Optics calculations do not take into account scattering effects due to 
sharp discontinuities in the geometry of an object. These discontinuities often account 
for a significant portion of the backscattered energy and should be included in the 
calculation. Different methods have been devised to deal with these discontinuities 
including the Geometrical Theory of Diffraction, the Uniform Theory of Diffraction, and 
the Physical Theory of Diffraction. The computer program used for the simulations here 
implements the Physical Theory of Diffraction, so that will be the focus of this 
discussion. Since the geometry files are created from triangles, every intersection of 
triangles creates an edge. Most of these intersections though are meant to approximate a 
smooth surface and not a diffracting edge. To ensure that scattering contributions are 
only calculated from diffracting edges, the angle between adjacent facets is checked and 
the edge is only used if it’s below a predetermined value. In the simulations here, that 
angle is chosen to be 120°. An illustration of scattering from an edge is given in Figure 8. 
A set of edges for the geometry is determined and again, only those edges which have 


incident wave illumination are used in the calculation. 


Source Angle: 165 





Wedge Angle: 30 


Figure 8. Scattering from a Diffracting Edge [19] 
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3. Shooting and Bouncing Rays 


The two previous methods only took into account the first bounce of the incident 
wave. After this initial bounce, the scattered waves can then hit other parts of the 
geometry and rescatter from them. If the electrical size of the target is large, then the 
bouncing of the waves can be approximated in the high frequency limit. Because the 
aircraft used in the simulation are much larger than the 3cm wavelengths used in the 
scattering, bouncing can be approximated using a method called Shooting and Bouncing 
Rays. Millions of rays spanning each of the different frequencies used are ‘shot’ at the 
target. As these rays encounter a facet of the target, they are reflected based on standard 
geometric techniques. The first facet encountered is ignored as the contribution from the 
first bounce is calculated with the previous methods. As the ray intersects another facet, 
the program checks to see if that facet’s scattered energy has a line of sight to the 
receiver. If so, the PO contribution from this facet based on the bouncing ray’s 
orientation is included in the scattered field. The ray then continues its bouncing path 
until the ray leaves the target geometry and no more bouncing can occur. Each facet that 
is encountered will contribute to the final scattered field. The limit on number of bounces 
used in the calculation was set at 20 bounces. In real targets, structures such as jet 
intakes can allow radar energy to bounce numerous times within the cavity with each 
bounce contributing to the cross section. Interior structures such as cockpits and engine 
compartments have limited model detail in the simulation so the detail provided from 
these structures is only marginally realistic. Interactions between wings, fuselage, 
weapons, and stabilizers are readily observed. A simplified visual representation of what 


the computer is actually doing is shown in Figure 9. 
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Figure 9. Rays Interacting With F-15 Model [18] 
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IV. CONTACT IDENTIFICATION 


Now that we know how the radar system determines that a contact has been made, 
the next step is determining what that contact is. The obvious answer of course is to just 
get on the radio and ask the pilot to identify himself. This method is included in the steps 
to identify a contact, but there are numerous problems associated with it. Problems 
include language differences, pilot deception, radio malfunction, operator error, and 
giving away too much information about the interrogator just by asking the question. 
The next most obvious answer is to visually look at the contact to determine what it is. 
This is a reliable method of identification but requires that friendly forces get within 
visual range of a contact which may or may not be friendly also. Visually identifying 
contacts becomes a dangerous and time-consuming process with the distances involved. 
Another rudimentary scheme is identification based on flight profiles. An aircraft flying 
standard commercial air routes with speeds and altitudes corresponding to a commercial 
airliner is most likely an airliner. This is of course an easily exploited method that the 
enemy can use for deception. Automated identification systems which can reach out to 
the range of the radar provide the best option for early warning of potential hostile 


contacts and many schemes have been devised and implemented over the years. 


A. IDENTIFICATION FRIEND OR FOE 


The most widely used and effective system is the Identification Friend or Foe 
(IFF) system [20]. It is an interrogation and reply based system with 2 channels of 
communication, one for the question and one for the reply (1.03 GHz and 1.09 GHz). 
Most air search radars are equipped with a separate antenna to send and receive IFF 
interrogations and replies. An example of one of these antennas is shown in front of the 


feed horn of the SPS-49 air search radar in Figure 10. 
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Figure 10. Example of IFF Antenna Attached to the SPS-49 Radar [21] 


Multiple modes of IFF are used, each with its own purpose. 


Mode 1 — The mission code used by military aircraft to identify their aircraft type 


and mission. 64 possible replies. 


Mode 2 — Used to identify the specific tail number of aircraft. There are 4096 


possible replies and only military aircraft use this mode. 


Mode 3/A — Mode used internationally for civilian air traffic control and identifies 
the aircraft’s flight number. It can also indicate special circumstances such as 


emergencies or hijackings. 


Mode 4 — Mode for encrypted military identification. Only with the correct 


cryptological gear and entered codes can this mode be successfully utilized. 


Mode C — Used for altitude broadcasting to provide air traffic controllers more 


information on the flight. 


IFF provides an efficient means for civilian air traffic controllers to direct planes 


in the air without having to fill the radio channels with requests for identification. It also 


gives military users the ability to determine if a contact is a friend. It does not give the 


ability to reliably distinguish a neutral contact like a civilian airliner from a hostile 


aircraft. To have effective military use, the current codes must always be entered into the 


hardware devices. Technician negligence or hardware malfunction make this system of 


identification less than perfectly reliable. Enemy eavesdropping of successful 
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interrogations may also allow them to reproduce the transmissions and masquerade as 
friendly aircraft although advancements in encryption technology significantly reduce 


this risk. 


B. NON-COOPERATIVE TARGET RECOGNITION 


When a friend or coworker comes in the room, we know immediately who that 
person is without having to ask them each time. Also when someone we do not know 
enters, but is wearing a U.S. Navy uniform, we know that person is probably not an 
enemy. We understand the identities of these people not by them telling us, but by 
recognizing features of their bodies, faces, and uniforms. Non-cooperative target 
recognition (NCTR) serves the same function of trying to identify a contact without 
asking who they are, and without them telling who they are. The goal is to have sensors 
like radar be able to tell what a contact is just by ‘looking’ at it. It is much more difficult 
for an enemy to change certain structural and mechanical characteristics of their aircraft 
to try to masquerade as something they are not. Various techniques have been developed 
to try to extract as much information as possible from the scattered radar energy in order 


to identify targets. 


1. Modulation Schemes 


Aircraft are dynamic machines in that they have constantly moving parts. Almost 
all aircraft have some sort of rotating machinery such as propellers, rotors, or fan blades. 
These rotating parts can interact with the incident radar pulse causing Doppler shifts at 
certain modulation frequencies based on the number, size, and rotation rate of the blades. 
When this principle is applied to aircraft, it is given the names Jet Engine Modulation 
(JEM), Propeller Rotor Modulation (PROM), and Helicopter Rotor Modulation (HERM) 
[22]. These signatures, shown as sidebands to the airframe return, can be compared to a 
library of known signatures to determine its aircraft type. For helicopters, this 
information is rather constant around the full azimuth range, but for planes, the frequency 
shifts change depending on azimuth and are most effective from the front or rear of the 


aircraft. This is a fairly effective method though if the modulation can be picked out of 
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the noise-laden radar return. The library required for this technique is not that extensive 


as only engine characteristics need to be stored. 


Z Inverse Synthetic Aperture Radar 


We often identify things based on pictures and Synthetic Aperture Radar (SAR) is 
an attempt to use radar waves rather than optical waves to create pictures [5]. Optical 
devices are able to achieve high resolution because the very small wavelengths used 


(microns) allow for a very small wavelength-to-aperture (1/A) ratio. The larger the 


aperture, the greater the possible cross range resolution. Since radar wavelengths are on 
the order of cm rather than microns, a much bigger aperture is required to achieve a 


similar 1/A ratio. By taking multiple discrete measurements while translating the radar, 


a larger effective aperture can be created. Inverse Synthetic Aperture Radar (ISAR) is 
mathematically related to SAR via a coordinate shift. The difference is: rather than a 
fixed object and a moving radar, the radar is stationary and the object is moving. 
Because of the wavelength size of radar waves, slight changes in range produce large 
changes of phase in the scattered field. Range must be measured very precisely and 
range fluctuations must be accurately accounted for. These pictures require some sort of 
relative rotation rate to acquire the necessary data so there is a time-lag associated with 
their generation. Once the pictures are created though, identification is as easy as looking 
at the picture. ISAR techniques are still being perfected and a sample image is shown in 


Figure 11. 
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Figure 11. Sample ISAR Image of a Boeing 727 airliner [5] 


3. Range Profiles 


The high range resolution techniques discussed in Chapter II enable identification 
methods such as range profiling [23]. Range profiles can be considered one-dimensional 
images as only down range data is presented with no cross range information. The 
profile is achieved by taking the inverse Fourier Transform of the frequency domain 
scattered field measurements. Converting the time domain into the range domain is a 
simple linear transformation r= ct. The greater the bandwidth of the receiver, the 


better the resolution in the range profile. 


A range profile shows peaks at the major scattering centers of the aircraft. The 
distance between peaks corresponds to the radial (range) distance between the centers. 
As the target angle of the aircraft changes, different parts of the aircraft will come in to 
illumination and change relative positions on the range profile. This means that a range 


profile can change significantly as the target rotates. Also, the multiple scattering 
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mechanisms that were discussed in Chapter III will cause peaks to show up in the range 
profile that do not correspond to a line of sight range, but rather the sum of the paths of 
the bounces. Also, parts of the aircraft will scatter strongly at some angles, and weakly at 
others. For instance, the wings of the aircraft may couple nicely with a horizontally 
polarized radar pulse where the vertical stabilizers will not. Conversely, if the aircraft 
banks 90 degrees, the stabilizers will be more effective scatterers and may show up 
strongly. A sample range profile of a model F/A-18 used later in the simulation is shown 


in Figure 12. 


The process of identification based on range profiles is a similar process to the 
modulation schemes previously discussed. A library of range profiles must be compared 
to the sampled profile, and which ever profile more closely matches is the likely identity. 
A library must include every aircraft expected to be encountered or the search will 
obviously report an incorrect match. Also, depending on the resolution of the system and 
how rotationally sensitive the profiles are, the library will have to have profiles covering 
fractions of degrees for every possible combination of azimuth and elevation. Another 
factor is that many aircraft are not configured the same for every flight. Some flights 
require external fuel tanks, bombs, missiles, etc. These extra features can significantly 
alter the range profile and they must also be taken into account in building the library. 
The libraries can be built with actual model aircraft on a range or with computer 
generated models and radar cross section estimation tools such as the ones used in the 
simulation. With all these factors to consider, obviously building dense enough libraries 
to reliably identify contacts will take millions of stored profiles. Searching the library 
must also be done in a reasonable amount of time (a few seconds) or the information 


won’t be actionable. 
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Figure 12. Simulated Range Profile of an F/A-18 


4. Bispectrum Profiles 


Range profiles contain scattering information for not just the first bounce, but for 
multiple bounces also. The information is depicted in the same manner though and it is 
almost impossible to visually separate the first bounce peaks from the multi-bounce 
peaks. Creating what are called bispectrum! images can help in identifying the multiple 
interactions of a target by separating them from the other terms [24]-[28]. The idea is 
that different aircraft types will have different characteristic interactions between the 
structures on the airframe. If these interactions can be isolated by the bispectrum, then a 


more informed comparison can be made with a library of signatures. 


The bispectrum is defined as the two-dimensional Fourier transform of the third 
order cumulant function of the scattered field measurements and is given by the 


following equation in the time domain where h(t) is the impulse response of the target: 


Bt, t,) = h(t, h(t, )h* (t, +t,) (12) 


! Also called birange or bitime, depending on the domain. I will use the term bispectrum to refer to 
this class of images. 
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The bispectrum therefore will have a strong peak wherever two peaks in the range profile 
can be associated with a peak further down in the profile. This indicates a possible 
interaction between these two scatterers. As each aircraft will have different interactions, 
these spots will show up at different places in the bispectrum, and a library can be 
searched in the same manner as before. It has been suggested that the bispectrum of a 
target has more rotational invariance than a range profile and therefore the library can be 
sparser. Bispectra also have the ability to suppress Gaussian white noise and therefore 
may be more useful in noisy environments than range profiles which do not suppress 


noise. 


Cy 3 POINT SCATTERERS EXAMPLE 


In order to illustrate the previous two methods, consider a simple widely used 
example — 3 point scatterers in space. It was presented in a previous thesis on the topic 


by Jiunn Wah Yeo [11]. The setup of the example is shown in Figure 13. 
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Figure 13. Illustration of 3 Point Scatterers Example [11] 


Yeo shows that when R,= 0m, R,=1.5m, and R,=3.6m, the scattered field can be 
reduced to 


an(1s)£ 4n(3.6)2 4a(s.)e| 


eeexen aga Plaae (13) 
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where C,,,,are constants. The Fourier transform of this signal will provide the range 


profile of these 3 point scatterers shown in Figure 14. The first 3 peaks correspond to the 
positions of the 3 points, but the re peak, noted as a “ghost” artifact, is due to the 


multiple scattering between the 3 points. 
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Figure 14. Range Profile of 3 Point Scatterers [11] 


A bispectrum plot can be created from the same frequency signal and is shown in 
Figure 15. In this image, the strong peaks indicate interactions by the second and third 
point. We also notice the symmetry properties of the bispectral image. A diagonal line 
through the origin in this particular image is the symmetry axis. A more complex target 
will naturally provide for more interactions and we will notice in the simulation, that 


much more complex bispectrum images are created by aircraft targets. 
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Characteristic Bispectrum Plot of 3 point Scatterers [11] 
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V. SIMULATION 


The goal of this thesis is to determine whether a classification system based on 
bispectrum images is a valuable enhancement or substitute for a system using range 
profiles. Since both systems need a library of aircraft signatures with which to compare, 
an experimental library had to be created. Afterwards, 4 sets of test cases were created 
and identifications of these test cases were made based on the two methods. A 


comparison of the simulation results follows in Chapter VI. 


A. BUILDING THE LIBRARY 
1. Scenario 


Since I’m not concerned with actually building a fieldable system and, instead, 
only with a performance comparison between algorithms, I do not need to have a library 
of every possible aircraft, configuration, and aspect angle. I begin with considering the 
situation I want to model. The most useful scenario for an early warning identification 
system is when a ship or aircraft detects an incoming contact at a distant range. Is the 
contact a potentially hostile enemy fighter, or is the contact a neutral commercial airliner? 
The sooner this question can be answered, the more likely we will be able to avoid 


another Vincennes-type disaster. This scenario is depicted in Figure 16. 


aw!) 


Figure 16. Simulation Scenario 
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2. Aircraft 


A real library of radar signatures should include every aircraft type that can be 
expected to be encountered on the mission. To build my test library, I included 10 
aircraft from a collection that would be similar enough to test the recognition methods, 
but also a few outside the class to provide a reality check. I chose eight fighter/attack 
aircraft, one helicopter, and one commercial airliner. The aircraft are from a collection of 
commercial 3-D models from a company called Mesh Factory [29]. The aircraft’s 
smooth surfaces are constructed from individual triangles and these triangles are what 
enable the use of the scattering techniques described in Chapter III. For instance, the 


F/A-18E Super Hornet, constructed from over 83,000 triangles, is depicted in Figure 17. 





Figure 17. Model and Wireframe of F/A-18E Super Hornet 


Fighter aircraft do not fly every mission with the same configuration. Some 
missions require external fuel tanks and an assortment of bombs, while others require 
only air-to-air missiles. In order to test the abilities of these methods to deal with 
alternate configurations, each of the fighter aircraft were also given a set of appropriate 
munitions for an additional test set, a ‘dirty’ test set. These additional configurations 


were not included in the library, only the test cases. The position of the helicopter’s main 


30 


rotor was also changed because every radar pulse on a helicopter will encounter the 
blades in a slightly different aspect allowing different scattering interactions. The Airbus 
A320 was not reconfigured as no external weapons seemed appropriate for the airframe. 


Appendix A provides a visual depiction of the aircraft and configurations used. 


a; Aspect Angles 


There is, of course, a performance/time trade-off between having a library dense 
enough to guarantee accurate identification, and a library sparse enough to be able to 
search through very quickly. In order to limit the size of the library, I first make the 
assumption that the aircraft are symmetric about the centerline — that is, that their port 
side is the same as their starboard side. This is not true of real aircraft as different 
munitions may be attached to each side and sensors may be attached anti-symmetrically. 
I neglect these effects for simplicity. A spherical coordinate system is used rather than 
the azimuth/elevation system. Because I am considering incoming contacts, I chose a 
range of target angles that would cover the expected vantage points. The angles are 
depicted in Figure 18. Each aircraft was scaled so that its wingspan is equivalent to the 


actual aircraft’s wingspan in meters. 
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Figure 18. Range of Angles Considered 


I gathered radar data for angles from 60° below vertical to 120° below vertical in 


5° increments, and from 0° ahead to 45° to port (315° relative) in 3° increments (which, in 


eal 


effect, covers from +45° to -45° due to symmetry). These angle combinations give 208 
profiles per aircraft providing a total library of 2080 searchable profiles. Admittedly this 
is a rather sparse library for extremely accurate identification, but it provides the means 


to compare the performance of the two methods side-by-side. 


4. Radars 


I consider radar contact acquisitions from aircraft radar and from shipboard radar 
to compare the differences in the performance of these methods at different frequencies 
and bandwidths. In order to ensure that this simulation remains unclassified and freely 
distributable, actual radar system specifications were not used. For the aircraft radar I 
simulate an X-band radar with 1 GHz bandwidth centered at 9.5 GHz. Using this 
bandwidth in equation (4) gives a range resolution of 15cm. For the shipboard radar I 
simulate a C-S band radar with a 400 MHz bandwidth centered at 4 GHz. This 
bandwidth provides a range resolution of 37.5cm. A separate library for each radar was 


built with the different scattering signatures. 


5. RCS Tools 


The radar cross section (RCS) estimation tool used to create the library is a 
program called Lucernhammer from Tripoint Industries [18]. It takes as input a facet file 
containing the aircraft geometries and calculates the scattered field using the PO, PTD, 
and SBR techniques described in Chapter III. The ability to model materials with 
different reflectivity and dielectric coefficients exists, but was not employed. The facet 
files are considered complete PEC objects which introduces inaccuracies in the 
signatures, but is satisfactory for comparison between the two methods. It generates a 
file containing the complex-valued scattered field measurements for each combination of 


@ and @ used for the library. The data are normalized so that the square of the 


magnitude is the radar cross section in square meters. 256 sampling points were taken 
within each bandwidth range allowing the range profiles to include the entire physical 
range of the aircraft models. The models were scaled to the appropriate dimensions 
mentioned before so that the wavelength relative to structure would be realistic. 


Lucernhammer has the ability to calculate scattered field in multiple polarizations, but the 
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horizontally transmitted and horizontally received polarization was chosen for the 
simulation. A special purpose file conversion utility was written to extract the desired 


information from the scattered field file and organize it into a Matlab friendly structure. 


6. Processing of Frequency Domain Data 


The frequency domain scattered field files were loaded into Matlab and separated 
into individual cell arrays. Each frequency array was normalized for unit energy. Range 
profiles were generated by applying the inverse Fourier transform to the simple 
rectangular windowed data. These range profiles cover a range of +/- 19.112 m and are 
256x1 arrays of complex-valued data. A library of bispectrum data was also generated 
using the code adapted from [11] in Appendix C. Each bispectrum is a 511x511 array of 


complex-valued data. 


B. RUNNING THE SIMULATION 
1. Test Cases 


The likelihood of trying to compare a test signal that is exactly the same as 
something in the library is practically non-existent. A radar return, especially from a 
distant target, is full of noise at different levels. Also, the aircraft will have slightly 
different configurations and will be a different aspect angles than what is in the library. I 
do not consider the effects of noise in the simulation. Since I was working with a rather 
sparse library, I was able to compare classifications at angles very near library angles and 
further away. My first data set was generated using the same aircraft configurations I 
used in building the library. I looked at combinations of angles 1 to 2 degrees away from 
library angles. I took data from each aircraft at 5 angle combinations spread across the 
library range and repeated those same combinations for each aircraft giving 50 total test 
cases for the first set. The second set is using the same aircraft configurations as the first, 
but with closer angle combinations within 1 degree of the library angles. 50 additional 
test cases were generated for this set. The purpose of these two sets is to compare the 
rotation invariance of the two identification methods and to see which method is 


degraded more by having a sparse library. 
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For the second two test case sets, the armed aircraft configurations described in 
Section V.A.2 were used. The first set includes exact angle combinations that are found 
in the library. Since noise was not added to the simulation, both methods would locate an 
exact replica of a library signal 100% of the time. Changing the configuration of the 
aircraft alone enables an analysis of the ability of the two methods to handle alternate 
configurations. 45 test cases for this set were created because the A320 was not 
reconfigured. The next set of test cases again use the armed configurations, but this time 
use the same close offset angles that were used in the second set of the stock 
configurations. In this set, the ability to handle both offset angles and multiple 
configurations is tested. Additionally, a direct comparison can be made with the second 
stock set to further determine how much configuration changes impact the identification 


methods. 


De Range Profile Correlation 


At the heart of the range profile based simulation is the question of how to 
determine which library entry most closely matches the test case. Because the profiles 
were all created from similarly normalized frequency data, they can be compared 
directly. The method of comparison used is the cross-correlation [5]. The correlation 


coefficient is the measure of how closely the two match and is given by the equation 
Cn, 0,9) = MAX | Py (Y Py + VO )dy" (14) 
where n is the type of aircraft and 0 and ¢ provide the aspect. The Cauchy-Schwartz 
inequality tells us that C(n,0,¢) <1. It equals one when the two signals are exactly alike 
(within a translation). The correlation giving the highest coefficient is selected as the 
identification and the likely aspect of the target is also extracted. The cross correlation is 
implemented in the simulation within the howCorrP.m file in Appendix C. It takes two 


complex range profiles and returns the correlation coefficient of the two. Matlab’s xcorr 


function is used with the magnitudes of the two signals. 
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a: Bispectrum Correlation 


The bispectrum based simulation also uses the idea of the cross-correlation but 
implements it differently [28]. Because the bispectrum profiles are 2-dimensional arrays, 
a 2 dimensional cross correlation must be used. The time or range domain based xcorr2 
function in Matlab would take an unreasonable amount of time to complete, so the 


correlation is done in the frequency domain. The equation used is 


HRA AIR A A) 
Y, 
ext] [ex 


fi fh 


Pino) = (15) 





Roh.) | 


R in this equation is the two-dimensional Fourier transform of the magnitude of the 
bispectrum profile. The correlation coefficient is extracted by taking the maximum value 
of this 2-D array over r, and r7. The implementation of this correlation can be found in 
the file howCorrBi.m in Appendix C. It takes as input the Fourier transform of two 
bispectrums and returns the correlation coefficient. Once again the identification is 


selected as the comparison with the highest correlation coefficient. 


4. Program Flow 


The correlation simulations operate by comparing each test case with every 
signature in the library and combining the results. The maximum overall correlation for 
each test case is designated the identification and noted in the output file. A matrix 
containing the tally for each identification versus test case group is also created. Test 
cases involving different aspects of the same aircraft are summed or averaged together. 
In addition to determining the maximum correlation for each test case, the simulation 
determines how closely each test case correlates to each aircraft type by filling an array 
with the maximum coefficient of correlating the test case with aircraft type. These 1x10 
arrays of coefficients for each test case are grouped again by test type and averaged to 
simulate the performance of averaging multiple radar hits together for increased 
identification accuracy. Simulation code is given in Appendix D and program flow logic 


flow is given below. 


35 


for each set of test cases 
load test case set 
open output file 
for each test case 
check test case with entire library 
make array out of all correlations 
find maximum correlation 
designate it the identification 
find maximum of subarrays 
add to coefficient matrix 
write identification to output file 
end for 
average coefficient matrix 
write matrices to output file 
end for 
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VI. ANALYSIS 


A. SIDE BY SIDE COMPARISON 


The similarities of the bispectrum plot and the range profile are evident when we 
look at them side by side with the same range scale as shown in Figure 19. The contour 
plot depiction of the bispectrum shows the highest peaks as red and lower points as blue. 
If we look at a horizontal line across the bispectrum at Range,=0, we see that the spikes 


in the range profile closely correspond to the peaks of the bispectrum. 
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Figure 19. Range Profile and Bispectrum Plot of A320 
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Comparing the range profile to the aircraft we see that the peak at a range of -12m 
likely corresponds to returns from the engine intake. If we look at a horizontal line (or 
vertical due to symmetry) at a range of -12m, we see many other peaks down the line. 
These peaks indicate that the scattered energy from the engine interacts with many 


structures down range — likewise with other structures. 


B. THE EFFECT OF BANDWIDTH 


As shown earlier, greater bandwidth allows greater range resolution. When trying 
to distinguish small fighter aircraft from one another, good range resolution is very 
important. The level of range profile detail lost from a front aspect view of an F-15 when 
we go from 1GHz bandwidth high frequency radar to the 400MHz bandwidth lower 
frequency radar is shown in Figure 20. The range resolution in the latter is not good 


enough to distinguish between parts of the aircraft very near others. 
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Figure 20. Range Profiles of F-1SE at HBW and LBW 
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Fewer features are present to enable us to distinguish the airplane from other 
similar fighters and we will see how this limitation will impact the classification 
programs. Just as the lack of bandwidth decreases the fidelity of the range profile, it also 
impacts the bispectrum plot. The difference in bispectrum plots of the same F-15 at the 
two bandwidths is shown in Figure 21. The tight peaks of the multiple interactions of the 
high bandwidth plot are much wider in the low bandwidth plot. The lower fidelity of the 
low bandwidth bispectrum plots will enable greater correlations with slightly different 


plots, and should cause classification errors. 
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Figure 21. Bispectrum Plots of F-15E at HBW and LBW 
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C. ROTATIONAL VARIANCE 


Slight rotations of the aircraft can significantly change the scattered field. The 
sharpness of the high bandwidth plots shows this effect very clearly. The difference in 
the range profile of an F/A-18E from @= 95° as we rotate from ¢= 21° to 27° is shown in 
Figure 22. We see how the many scattering points of the engines and weapon pylons 
interact and change with a small rotation. This sensitivity highlights the importance of 


having a library with very close angles and the simulation will further illustrate this. 
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Figure 22. Rotational Variance of F/A-18E Range Profiles 
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High Bandwidth Bispectrum of F/A-16E at Theta=95 Phi=21 
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Figure 23. Rotational Variance of F/A-18E Bispectrum Plots 


The same F/A-18 aspect and rotation is shown in Figure 23. We see that even 
with the bispectrum, the plots change significantly. The radar wave bouncing from 
structure to structure is very sensitive to incident angle. By the time the aircraft rotates 
the full 6°, the plot looks very dissimilar to the original. Tables 1-4 show this effect with 
correlation coefficients. The model F/A-18E at 9= 95° and ¢= 21° was correlated with 


its neighbors using range profiles and bispectrums at different bandwidths. The profile 


4] 


correlations at each bandwidth returned similar results and we note that the rotation does 
not change the range profile correlation by a large percentage. This may cause problems 
when trying to classify similar targets. The bispectrum correlations show that the test 
case did not correlate well with its neighbors. A large percentage change in the 
correlation coefficient is returned when the contact is rotated by a few degrees. We also 
note that the low bandwidth provided less of a difference than the high bandwidth as 


expected. This large difference in rotational similarity should not make the bispectrum 


able to handle large angle disparity. 
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Table 1. | Rotational High Bandwidth Profile Correlations 
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Table 2. Rotational Low Bandwidth Profile Correlations 


Phi 
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95 | 1.0000 | 0.0859 | 0.1236 
100} 0.0991 | 0.0783 
105|{ 0.0880 | 0.0940 | 0.0934 


Table 3. Rotational High Bandwidth Bispectrum Correlations 
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Table 4. Rotational Low Bandwidth Bispectrum Correlations 
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D. CONFIGURATION CHANGES 


A change in the configuration of an aircraft has a great effect on the range profile 
and bispectrum. External fuel tanks, missiles, and bombs can contribute a large portion 
of the backscattered signal. The difference in the range profile of an unarmed F-16 with 
O= 120° and ¢= 45° and the armed version at the same aspect is shown in Figure 24. 
The additional structures cause peaks to appear in the profile where none exist in the 
unarmed profile. The bispectrum plots of the same aircraft are shown in Figure 25. We 
see that not only do new scattering interactions appear, but interactions present in the 
unarmed plot disappear. This can be due to masking of structures by the armaments or 


interference of the free path of the waves. 
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High Bandwidth Range Profile of Armed F-16C at Theta=120 Phi=45 
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Figure 24. Multi-Configuration Range Profiles of F-16C 
43 


High Bandwidth Bispectrum of Unarmed F-16C at Theta=120 Phi=45 


Range (m) 





Range (m)} 


High Bandwidth Bispectrum of Armed F-16C at Theta=120 Phi=45 
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Figure 25. Multi-Configuration Bispectrum Plots of F-16C 


E. SIMULATION RESULTS 

The results of the simulations described in Chapter V are analyzed by section 
below. The individual results are provided in Appendix E. 

1. High Bandwidth Results 


The high bandwidth (HBW) range profile simulations involving the stock aircraft 


showed promising results. Based on individual test cases, the success rate was only 60% 
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per set, but when multiple radar hits of the same test aircraft were averaged together, the 
highest correlation coefficient corresponded to the correct aircraft every time. The 
highlighted block in the tables indicate the maximum coefficient in each row. There is 
not much margin for error though as the coefficients all range between 0.75 and 0.90. 
The hypothesis before the tests were run was that the second set, with angles very near 
the library angles, would show a greater correlation than the first set. Tables 5 and 6 
show that hypothesis was incorrect for the range profiles — it was true for only 5 of the 


aircraft and the coefficients had a similar range in each table. 
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Table 5. |HBW Range Profile Simulation With Stock Aircraft at Far Angles 


Library Aircraft 


| CF AS20 [EAGB] F4 | Fi4 | F15 | F16 | FAI8| MiG _| Su27 | UH6O | 
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Table 6. HBW Range Profile Simulation With Stock Aircraft at Near Angles 
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The same two test sets run with the bispectrum algorithm produce different 
results. There was a 58% individual success rate for the first set and an 80% success rate 
for the second set with closer angles. Because the HBW bispectrum plots produce such 
narrow peaks, the plots are more easily misaligned and the correlation coefficients are in 
the .4-.7 range. What is important is the difference between the maximum coefficient 
and the others in the row. After averaging the test cases, the simulation misidentified the 
A320 as an EA6B in the first set of tests. The rest of the test cases were correctly 
identified after averaging. Because of the greater individual success rate and the greater 
difference between coefficients in a row, the bispectrum at near angles outperforms the 


range profiles with stock aircraft. 


Library Aircraft 


| | A20 | EASB] F4 | Fi4 | F15 | F16 | FAI8 | MiG | Su27 | UH6O | 
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Table 7. | HBW Bispectrum Simulation With Stock Aircraft at Far Angles 
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Table 8. | HBW Bispectrum Simulation With Stock Aircraft at Near Angles 
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Now let’s look at how the two methods performed with the armed aircraft test 
sets. The first test set includes armed aircraft at angles included in the library and the 
second set used the same angles as the second stock set. Tables 9 and 10 show the 
results. The first row includes zeroes because there were no A320 test cases. The range 
profile simulation with the first set produced an individual success rate of 73%. After 
averaging, the simulation chose the correct identification every time with the first set with 
high correlation coefficients. The EA-6B, MiG-29, and Su-27 had the highest 
coefficients of the set as those aircraft were armed with only a couple extra missiles. The 
other aircraft were changed more with the addition of a variety of ordinance and external 
fuel tanks, and their coefficients were reduced as a result. The second set produced an 
individual success rate of only 27% and after averaging, only 4 of the aircraft were 
correctly identified. These aspects had less than 1° of deviation from the library but the 


extra change provided by the configurations caused more stress on the simulation. 


Library Aircraft 


| | A20 | EAGB] F4 | Fi4 | F15 | F16 | FAI8 | MiG | Su27 | UH6O | 
0. 3365 


Test Cases 


; 8167 OBIT osora|o ane 03139 SEE 0.9405] 0. 
UH60 |0.8252|0.8102| 0.8715] 0.8303 9] 0.8652] 0.8506 | 0.8359| 0.8306] 0.879 





Table 9. HBW Range Profile Simulation With Armed Aircraft at Exact Angles 
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Table 10. HBW Range Profile Simulation With Armed Aircraft at Near Angles 
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The bispectrum simulation was able to handle these two test sets more easily. 
The first set produced an individual success rate of 80%. After averaging, the simulation 
correctly identified every aircraft as did the range profile simulation. The effect of the 
varying levels of configuration change is more easily noticed with the coefficients 
produced by the bispectrum. The slightly changed aircraft maintained coefficients 
around 0.9 contrasting with the most heavily armed F/A-18 with a coefficient of only 
0.57. The second test set was again harder to manage for the simulation. It produced an 
individual success rate of 49%, but much greater than the 27% for the range profile based 
target identification. After averaging though, 7 of the 9 aircraft were correctly identified, 
thus outperforming the range profile’s 4 identifications. In the two test rows that were 
misidentified, the correct aircraft came in second place by a small margin. Tables 11 and 


12 contain the results. 


Library Aircraft 
| | AS20 | EASB] F4 | Fi4 | F15 | F16 | FA18| MiG | Su27 | UH6O | 
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Table 11. . HBW Bispectrum Simulation With Armed Aircraft at Exact Angles 
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Library Aircraft 
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Table 12. .HBW Bispectrum Simulation With Armed Aircraft at Near Angles 


2s Low Bandwidth Results 


The reduction in the ability to distinguish structures on the aircraft caused the low 
bandwidth (LBW) simulations to produce worse results. The range profile simulations 
with the stock aircraft configurations produced an individual success rate of 40% for the 
far angles and 48% for the near angles. Both these are lower than the 60% for the HBW 
cases. After averaging, the simulation failed to correctly identify the F-15 for the far 
angles. Curiously, with the near angles, it failed with the F-14, F-16, and MiG-29. The 
difference must lie in the particular angles chosen for the different test sets. Since the 
LBW cases are not as rotationally variant, there should be less difference between the 
near and far angle cases. Despite having a better individual success rate, the second set 
underperformed the first set after averaging. The troublesome issue with these 
simulations was the very narrow margin between the correct and incorrect identifications. 


Tables 13 and 14 show the results. 
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Table 13. LBW Range Profile Simulation With Stock Aircraft at Far Angles 
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Table 14. |. LBW Range Profile Simulation With Stock Aircraft at Near Angles 


The bispectrum simulation with the first test set produced an individual success 
rate of 56% outperforming the range profile, but still failed to identify the F-15 after 
averaging. The success rate is on par with the 58% achieved with the HBW set. The 
second test set with near angles produced an individual success rate of 42% while 
misidentifying two aircraft after averaging. As with the range profile simulation, the 
bispectrum simulation was not as effective with near angles as with far angles at LBW. 
The LBW did produce greater correlation coefficients than with the same test cases at 
HBW. Tables 15 and 16 show the results. As expected though, the HBW stock cases for 
both the bispectrum and the range profile outperformed the LBW simulations. 
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Table 15. | LBW Bispectrum Simulation With Stock Aircraft at Far Angles 
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Table 16. LBW Bispectrum Simulation With Stock Aircraft at Near Angles 





The hypothesis was that configuration changes would have less effect with LBW 
cases because of the reduced ability to pick out the changes in the profiles. This 
generally held true. The range profile simulation with the first test set of armed aircraft 
produced an individual success rate of 71% and misidentified one aircraft after averaging. 
Correlation coefficients were all higher compared to the same test set at HBW. The 
second test set produced a success rate of 36% and misidentified 3 of the 9 aircraft. This 
36% was only slightly worse than the second stock set with the same angles. The drastic 
difference of the same two test sets at HBW shows that configuration changes have less 
effect at LBW. Whereas the HBW range profiles outperformed the LBW sets with stock 
configurations, the second LBW armed set outperformed the same HBW set. Tables 17 


and 18 show the results. 
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Table 17. | LBW Range Profile Simulation With Armed Aircraft at Exact Angles 
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Table 18. LBW Range Profile Simulation With Armed Aircraft at Near Angles 


The bispectrum simulations with armed aircraft followed closely with the same 
range profile simulations. The first test set had a success rate of 73% and misidentified 
one aircraft after averaging. This bispectrum set produced a bigger difference between 
the coefficients for the correct and incorrect aircraft than the same range profile set. This 
creates a stronger likely identification. The second test set had a success rate of 40% and 
misidentified 2 aircraft slightly outperforming the same range profile set. This second 
armed test set closely matches the performance of the second stock test set again showing 
that the LBW simulations are less affected by configuration changes. The LBW 
bispectrum armed test sets only perform slightly worse than the same sets at HBW. 


Tables 19 and 20 show the results. 
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LBW Bispectrum Simulation With Armed Aircraft at Exact Angles 
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LBW Bispectrum Simulation With Armed Aircraft at Near Angles 
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Vil. CONCLUSIONS 


A. SUMMARY 


Military members make life or death decisions based on the set of information 
they have at their disposal at the instant the decision needs to be made. A basic 
component of that information is the question of whether the contact about to be 
destroyed is an enemy or a friend. Any hesitation made by the operator as a result of not 
trusting the identification being presented by his system can be fatal. The goal is to 
maximize the quality of the information he receives in order to build the trust that when 
the system tells the operator the contact is an incoming F-16, he knows that it is. In 
instances where this information was not present, or where it was incorrect, the resulting 


consequences were shown. 


A radar system detects contacts by sending and receiving electromagnetic energy 
and processing the received signal. The basic hardware layout was described and how 
radar parameters such as beam width, pulse width, bandwidth, and frequency affect the 
ability of the system to resolve targets. The electromagnetic energy interacts with the 
target and scatters to form unique signatures. The RCS codes used in the simulation 
estimate the scattered field by summing contributions from physical optics, the physical 
theory of diffraction, and multiple interactions. By analyzing the scattered field we can 
try to determine what kind of contact created the field. A library of various aircraft was 
created to facilitate the comparison of the range profile and bispectrum methods of 


identification. 


B. EFFECTIVENESS 


It was shown through a series of simulations that a system using bispectrum 
correlation meets or exceeds the performance of a system using range profile correlation. 
The bispectrum produced a 5.3% increase in identification accuracy over range profiles 
after averaging in both the HBW and LBW simulations. The results are summarized in 
Tables 21 and 22. Future sensor systems should have an increased bandwidth over 
current systems which make techniques such as the bispectrum more promising. 
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| [Range Prof 


Individual 55.0% 66.8 % 
86.8 % 92.1 % 


Table 21. Results Summary for High Bandwidth Simulation 


| [Range Prof 





Individual 48.8% 52.8% 
78.9 % 84.2 % 


Table 22. Results Summary for Low Bandwidth Simulation 





The various success rates shown throughout the thesis are not impressive, but they 
should not be considered as stand-alone figures. The library used to search through was 
very coarse and is only useful to have a standard to compare the two algorithms side-by- 
side. Also, much more advanced signal processing techniques can be used to refine the 
methods. Standard methods, initial conditions, and test cases were employed in the 


simulation, allowing direct comparison between the methods. 


No one method should be used for identification; rather, many methods can be 
combined in a hybrid system. For instance, since the range profile system takes 
significantly less time to run, the field of possible matches can be quickly narrowed using 
range profiles, and then the likely matches can be further analyzed with the bispectrum. 
Also, weighting factors can be applied based on the likelihood of that contact actually 
being present in a certain region. Many clever techniques have been developed for 
identification, but this thesis shows that using the bispectrum should remain a part of that 


set. 


c. ADDITIONAL CONSIDERATIONS 


More factors go into the determination of a good identification method than were 
considered in this simulation. Implementation issues such as speed of execution, 
efficiency of the algorithms, and size of storage required to hold the library should be 
also considered. For instance, the bispectrum simulation used in the thesis took over six 
minutes per test case to search through the entire library on a modern computer. The 


range profile simulation took just over six seconds per test case. Additionally, the library 
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of range profiles took up only a couple of megabytes of storage space whereas the library 
of bispectra took up three orders of magnitude greater than that. Correlating with a data 


set that would have to be much larger than this is not a trivial process. 


The effect of noise was also not considered. The very faint radar returns of 
distant targets are full of noise and can drastically affect the identification ability. It has 
been shown that the bispectrum suppresses additive Gaussian noise, which may give it an 
additional advantage over the range profiles [24]-[28]. The addition of different levels of 
noise into the test cases should be considered to determine the signal-to-noise ratio when 


these methods break down. 


The polarizations used for the returns in this thesis were the horizontally 
transmitted, horizontally received polarizations. Vertical or cross polarizations should be 
considered to see if the bispectrum is more useful with multiple polarization information. 
It may be that multiple interaction scattering tends to produces certain polarizations and 


that could be exploited with the bispectrum. 


Clutter was added into the simulation in the form of configuration changes. A 
real library might have multiple standard configurations of aircraft already considered so 
additional clutter should also be added. Whether this method would be useful in clutter 
filled areas such as against vehicles on the ground or ships on choppy waters should be 
determined. The current styles of ships in the fleets of the world have so many external 
structures to facilitate multiple interactions that the bispectrum would likely be a useful 


tool to identify distant ships. 
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APPENDIX 


A. AIRCRAFT USED IN SIMULATION 





Figure 26. Airbus A320 
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Figure 27. EA-6B Prowler 





Figure 28. F-4N Phantom II 


er - 


Figure 29. F-14D Tomcat 


aa 


Figure 30. F-15E Strike Eagle 








Figure 31. F-16C Falcon 





Figure 32. F/A-18E Super Hornet 
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Figure 33. MiG-29A Fulcrum 





Figure 34. Su-27 Flanker B 





Figure 35. UH-60L Black Hawk 
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B. 


LIST OF TEST CASES 


a320 
a320 
a320 
a320 
a320 
ea6b 
ea6b 
ea6b 
ea6b 
ea6b 
{4 
f4 
f4 
f4 
4 
f14 
f14 
f14 
f14 
f14 
f15 
f15 
f15 
f15 
f15 
f16 
f16 
f16 
f16 
f16 
fal8 
fal8 
fal8 
fal8 
fal8 
mig29 
mig29 
mig29 
mig29 
mig29 
su27 
su27 
su27 
su27 
su27 
uh60 
uh60 
uh60 
uh60 
uh60 
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Dirty2 
8 ® 
100 =000.5 
081 009 
086 = =024.5 
100 = 031 
091 040 
100 000.5 
081 009 
086 = 024.5 
100 = 031 
091 040 
100 000.5 
081 009 
086 = 024.5 
100 = 031 
091 040 
100 = 000.5 
081 009 
086 = =024.5 
100 = 031 
091 040 
100 = 000.5 
081 009 
086 =024.5 
100 = 031 
091 (040 
100 = 000.5 
081 009 
086 = =024.5 
100s 031 
091 040 
100 = 000.5 
081 009 
086 = =024.5 
100 = 031 
091-040 
100 = 000.5 
081 009 
086 = 024.5 
100 = 031 
091 040 
100 = 000.5 
081 009 
086 = =024.5 
100 = 031 
091-040 


c. BISPECTRUM CODE 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% This function calculates and returns the bispectrum of the input signal. 
% function adapted from code by Jiunn Wah Yeo [11] 


% signal is the freq domain complex scattered field 


function output = getBi(signal); 


% generate third order cumulant of signal 
signal_length = length(signal)-1; 


% Assignment and pre-conditioning of signal 
signal_conj=conj(signal(:)); 
signal_col=signal(:); 


N = length(signal_conj); 


% Triple correlation of signal 
X1 = repmat(signal_conj,1,2*signal_length+1); 


C2=[zeros(signal_length,1);signal_col(1:signal_length+1)]; 
R2=[signal_col(signal_length+1:N).' zeros(1,signal_length)]; 
X2=hankel(C2,R2); 

X3=X2.'; 


cum3=X2*(X1.*X3)/N; % 3rd order cumulant 


% Perform Fourier Transform of the third order cumulant to 
% generate the bispectrum 


output = fftshift(fft2(ifftshift(cum3))); % bispectrum 


D. SIMULATION CODE 


% AS 2A he 2g 2 he 2g 2s fe 2g 2s 2g 2g 2 2k 2g oo 2 p_corr_sim.m BRS 2g 2k os 2g 2k os 2g 2k 2s 2k 2k os os 2k ois os 2k oi oi 2 ok 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% RUN TIME: 5 minutes 


% This file implements a cross correlation classification scheme 
% for range profiles. 


% This is the general program control flow for running the correlation 

% simulation. It is broken up into 4 main loops. Each loop cycles through 
% the a set of test cases based on 10 aircraft for stock configurations 

% and 9 aircraft for armed configurations according to the table below. It 
% creates 6 output files, one for each loop and 2 for consolidating the 

% different bin matrices. 


% inputs: 
% p_stock same configs as library at angles far from library (1-3 de 

ry g 
% p_stock2 same configs as library at angles closer to library (<=1 deg) 
% p_dirty armed aircraft at angles corresponding to library (0 deg) 
% p_dirty2 armed aircraft at same angles as stock2 
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% lookup_table a text array for displaying user friendly results 
% ../zLibrary/ holds the library of aircraft 


% class_array -- an array of 2080 correlation coefficients from comparing 
% one test case to the library. 
% index -- the maximum corr coef in the array 


% helper files: checkLibraryCORRP.m, rowS.m, rowD.m, column.m, 
% fillRowWithCoefs.m, howCorrP.m, ../zLibrary/ 


clear; 
load lookup_table; 


% these are the matrices for individual identifications 
bin_matrix_s = zeros(10,10); 

bin_matrix_s2 = zeros(10,10); 

bin_matrix_d = zeros(10,10); 

bin_matrix_d2 = zeros(10,10); 


% these are the matrices for averaged correlation coefficients 
coef_matrix_s = zeros(10,10); 

coef_matrix_s2 = zeros(10,10); 

coef_matrix_d = zeros(10,10); 

coef_matrix_d2 = zeros(10,10); 


out = fopen(‘p_corr_stock_results.txt’,'wt'); Yoopen file for results 

fprintf(out, ‘p_stock corr tests\n'); 

load p_stock; “load test case set 

for testcase = 1:50 
% check one test case with library 
class_array = checkLibraryCORRP(p_stock {testcase} ); 
% find most likely aircraft 
index=find(class_array==max(class_array)); 
% fill matrix with coefficients 
coef_matrix_s(rowS(testcase),:) = coef_matrix_s(rowS(testcase),:) + fillRowWithCoefs(class_array); 
% fill binary matrix with likely identification 
bin_matrix_s(rowS(testcase),column(index))=bin_matrix_s(rowS(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear p_ stock; 

coef_matrix_s = coef_matrix_s / 5; Y%average matrix (5 cases per aircraft) 

felose(out); 





out = fopen(‘p_corr_stock2_results.txt','wt'); Yoopen file for results 

fprintf(out, ‘p_stock2 corr tests\n'); 

load p_stock2; “load test case set 

for testcase = 1:50 
% check one test case with library 
class_array = checkLibraryCORRP(p_stock2 {testcase} ); 
% find most likely aircraft 
index=find(class_array==max(class_array)); 
% fill matrix with coefficients 
coef_matrix_s2(rowS(testcase),:) = coef_matrix_s2(rowS(testcase),:) + fillRowWithCoefs(class_array); 
% fill binary matrix with likely identification 
bin_matrix_s2(rowS(testcase),column(index))=bin_matrix_s2(rowS(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear p_stock2; 

coef_matrix_s2 = coef_matrix_s2 / 5; Y%average matrix (5 cases per aircraft) 

fclose(out); 





out = fopen(‘p_corr_dirty_results.txt’,'wt'); “open file for results 
fprintf(out, 'p_dirty corr tests\n'); 
load p_dirty; Yoload test case set 
for testcase = 1:45 
% check one test case with library 
class_array = checkLibraryCORRP(p_dirty {testcase} ); 
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% find most likely aircraft 
index=find(class_array==max(class_array)); 
% fill matrix with coefficients 
coef_matrix_d(rowD(testcase),:) = coef_matrix_d(rowD(testcase),:) + fillRowWithCoefs(class_array); 
% fill binary matrix with likely identification 
bin_matrix_d(rowD(testcase),column(index))=bin_matrix_d(rowD(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear p_dirty; 

coef_matrix_d=coef_matrix_d/5; %average matrix (5 cases per aircraft) 

fclose(out); 





out = fopen(‘p_corr_dirty2_results.txt','wt'); Yoopen file for results 

fprintf(out, ‘p_dirty2 corr tests\n'); 

load p_dirty2; “load test case set 

for testcase = 1:45 
% check one test case with library 
class_array = checkLibraryCORRP(p_dirty2 {testcase} ); 
% find most likely aircraft 
index=find(class_array==max(class_array)); 
% fill matrix with coefficients 
coef_matrix_d2(rowD(testcase),:) = coef_matrix_d2(rowD(testcase),:) + fillRowWithCoefs(class_array); 
% fill binary matrix with likely identification 
bin_matrix_d2(rowD(testcase),column(index))=bin_matrix_d2(rowD(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear p_dirty2; 

coef_matrix_d2 = coef_matrix_d2 / 5; %average matrix (5 cases per aircraft) 

fclose(out); 





% Consolidate and print results to a file 

% results will be matrices where rows correspond to test case platforms (5 

% each) and columns correspond to library identifications by aircraft type. 

out = fopen(‘p_matrix_bin.txt’,'wt'); 

fprintf(out, ‘p_stock corr tests\n\n'); 

fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 

fprintf(out,'%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f \n', bin_matrix_s'); 
fprintf(out, '\np_stock2 corr tests\n\n'); 
fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 

fprintf(out, '%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f \n', bin_matrix_s2'); 
fprintf(out, '\np_dirty corr tests\n\n'); 
fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 

fprintf(out, '%1.0f %1.0f %1.0f %1.0f %L.Of %LOf %LOf %L.Of %1.0f %1.0f \n', bin_matrix_d'); 
fprintf(out, '\np_dirty2 corr tests\n\n'); 
fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 

fprintf(out,'%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f \n’, bin_matrix_d2'); 
fclose(out); 





% Consolidate and print results to a file 

out = fopen('p_matrix_coef.txt','wt'); 

fprintf(out, ‘p_stock corr tests\n\n'); 

fprintf(out,'A320 EA6B F4 F114 FI5 F16 FAI8 MiG Su27 UH60\n'); 

fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n', coef_matrix_s'); 
fprintf(out, '\np_stock2 corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F14 FI5 F16 FAI8 MiG Su27 UH60\n'); 

fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n', coef_matrix_s2'); 
fprintf(out, '\np_dirty corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F14 FI5 F116 FAI8 MiG Su27 UH60\n'); 

fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n', coef_matrix_d'); 
fprintf(out, '\np_dirty2 corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F14 FI15 F116 FAI18 MiG Su27 UH60\n'); 

fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n', coef_matrix_d2'); 
fclose(out); 

















% AS 2S AR 2S is 2g 2S 2s 2g 2 os 2g 2g os ok 2k oi ok 2 checkLibraryCORRP.m OAS sg 26 2S is 2g 2k 2s fe 2k 2s 2k 2k os ok 2k ok 2 2k 2 oi 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 
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% checkLibrary is a function that provides an array where each element 
% indicates how similiar the signal is to the specific library signal. 

% The array element corresponds to a lookup table of all the library 

% signals. 


function output_array = checkLibraryCORRP(signal); 


load ../zLibrary/a320_pc; “%load a320 section of library 

for index=1:208 for the a320 section of the output_array 
output_array(index) = howCorrP(signal,a320_pc {index}); 

end 

clear a320_pc; 


load ../zLibrary/ea6b_pc; Yoload ea6b section of library 

for index=209:416 %for the ea6b section of the output_array 
output_array(index) = howCorrP(signal,ea6b_pc {index-208}); 

end 

clear ea6b_pc 


load ../zLibrary/f4_pc; “load f4 section of library 

for index=417:624 %%for the f4 section of the output_array 
output_array(index) = howCorrP(signal,f4_pc {index-416}); 

end 

clear f4 pc 


load ../zLibrary/f14_pc; %load f14 section of library 

for index=625:832 %%for the f14 section of the output_array 
output_array(index) = howCorrP(signal,f14_pc {index-624}); 

end 

clear f14_pc 


load ../zLibrary/f15_pc; Y%load f15 section of library 

for index=833:1040 %for the f15 section of the output_array 
output_array(index) = howCorrP(signal,f15_pc {index-832}); 

end 

clear f15_pc 


load ../zLibrary/f16_pc; %load f16 section of library 

for index=1041:1248 %for the f16 section of the output_array 
output_array(index) = howCorrP(signal,f16_pc {index-1040}); 

end 

clear f16_pc 


load ../zLibrary/fal8_ pc; Yoload fal8 section of library 

for index=1249:1456 %for the fal8 section of the output_array 
output_array(index) = howCorrP(signal,fal8_ pc {index-1248}); 

end 

clear fal8_ pc 


load ../zLibrary/mig29 pc; “%load mig29 section of library 

for index=1457:1664 %for the mig29 section of the output_array 
output_array(index) = howCorrP(signal,mig29_pc {index-1456}); 

end 

clear mig29_pc 


load ../zLibrary/su27_ pc; Yoload su27 section of library 

for index=1665:1872 %for the su27 section of the output_array 
output_array(index) = howCorrP(signal,su27_pc {index-1664}); 

end 

clear su27_pce 


load ../zLibrary/uh60_pe; “%load uh60 section of library 

for index=1873:2080 %for the uh60 section of the output_array 
output_array(index) = howCorrP(signal,uh60_pc {index-1872}); 

end 

clear uh60_ pce 


66 


% AR 2g os fe 2A 2s fe 24g 2s ig 2g 2 2s 2k 2 oi 2k 2 ok 2 howCorrP.m ORG 8 2g 2S 2S is 2g 2S 2 2g 26 2S ok 2k 2 ok 2k 2g 2 2 2 OK 2 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% This file includes the definition of the howNear function. 
% howNear determines how closely the test signal matches the library 
% signal. 


% input: test signal <256x1> and library signal <256x1> 
% output: corr coef between 0 and | 


function result = howCorrP(testsig,librarysig); 


result = max(xcorr(abs(testsig),abs(librarysig)))*256; 


% AS 2 os hg 2s 2s 2g 2g 2s fg 2g 2s is 26 2 2k 2g oo 2 bi_corr_sim.m OAS sg AS 2S ig 28 2s ig 26 2s 2g 2S 2s 2k 2S os 2k 2 oR ok 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% RUN TIME: 40 Hours 


% This file implements a cross correlation classification scheme 
% for Bispectrum data. 


% This is the general program control flow for running the bispectrum 

% simulation. It is broken up into 4 main loops. Before each loop the 

% actual checking of the library is done with another function. Inside the 
% loop, the arrays produced are filtered for results. This control 

% organization differs from p_corr_sim in that the time complexity is 

% reduced by decreasing inefficient code. Time becomes an issue for this 
% code and expect this file to take 40 hours to run completely on lab 

% computers. This can be run in quarters without difficulty. It 

% creates 6 output files, one for each loop and 2 for consolidating the 

% different bin matrices. 


% inputs: 

% R_stock same configs as library at angles far from library (1-2 deg) 
% R_stock2 same configs as library at angles closer to library (<=1 deg) 
% R_dirty armed aircraft at angles corresponding to library (0 deg) 

% R_dirty2 armed aircraft at same angles as stock2 

% lookup _table a text array for displaying user friendly results 

% ../zLibrary/ holds the library of aircraft 


% class_array -- a cell array of arrays of 2080 correlation coefficients 
% from comparing each test case to the library. 
% index -- the maximum corr coef in the array 


% helper files: checkLibraryCORRBi.m, rowS.m, rowD.m, column.m, 
% fillRowWithCoefs.m, howCorrBi.m 


clear; 
load lookup_table; 


% these are the matrices for individual identifications 
bin_matrix_s = zeros(10,10); 

bin_matrix_s2 = zeros(10,10); 

bin_matrix_d = zeros(10,10); 

bin_matrix_d2 = zeros(10,10); 


% these are the matrices for averaged correlation coefficients 
coef_matrix_s = zeros(10,10); 

coef_matrix_s2 = zeros(10,10); 

coef_matrix_d = zeros(10,10); 

coef_matrix_d2 = zeros(10,10); 
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load R_ stock; “load set of test cases 

class_array = checkLibraryCORRBi(R_stock,50); Yoperform all checks 

out = fopen(‘bi_corr_stock_results.txt','wt'); Yoopen file for results 

fprintf(out, 'bi_stock corr tests\n'); 

for testcase = 1:50 
% find most likely aircraft 
index=find(class_array {testcase} ==max(class_array {testcase} )); 
% fill matrix with coefficients 
coef_matrix_s(rowS(testcase),:) = coef_matrix_s(rowS(testcase),:) + fillRowWithCoefs(class_array {testcase} ); 
% fill binary matrix with likely identification 
bin_matrix_s(rowS(testcase),column(index))=bin_matrix_s(rowS(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear R_ stock; 

coef_matrix_s = coef_matrix_s / 5; Y%average matrix (5 cases per aircraft) 

fclose(out); 





load R_stock2; “load set of test cases 

class_array = checkLibraryCORRBi(R_stock2,50); “operform all checks 

out = fopen(‘bi_corr_stock2_results.txt',wt'); Yoopen file for results 

fprintf(out, "bi_stock2 corr tests\n'); 

%sort array for results 

for testcase = 1:50 
% find most likely aircraft 
index=find(class_array {testcase}==max(class_array {testcase} )); 
% fill matrix with coefficients 
coef_matrix_s2(rowS(testcase),:) = coef_matrix_s2(rowS(testcase),:) + fillRowWithCoefs(class_array {testcase} ); 
% fill binary matrix with likely identification 
bin_matrix_s2(rowS(testcase),column(index))=bin_matrix_s2(rowS(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear R_stock2; 

coef_matrix_s2 = coef_matrix_s2 / 5; Y%average matrix (5 cases per aircraft) 

fclose(out); 





load R_dirty; Yoload set of test cases 

class_array = checkLibraryCORRBi(R_dirty,45); Yoperform all checks 

out = fopen(‘bi_corr_dirty_results.txt','wt'); Yoopen file for results 

fprintf(out, "bi_dirty corr tests\n'); 

for testcase = 1:45 
% find most likely aircraft 
index=find(class_array {testcase} ==max(class_array {testcase} )); 
% fill matrix with coefficients 
coef_matrix_d(rowD(testcase),:) = coef_matrix_d(rowD(testcase),:) + fillRowWithCoefs(class_array {testcase} ); 
% fill binary matrix with likely identification 
bin_matrix_d(rowD(testcase),column(index))=bin_matrix_d(rowD(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear R_dirty; 

coef_matrix_d=coef_matrix_d/5; Y%average matrix (5 cases per aircraft) 

fclose(out); 





load R_dirty2; “load set of test cases 

class_array = checkLibraryCORRBi(R_dirty2,45); Yoperform all checks 

out = fopen(‘bi_corr_dirty2_results.txt','wt'); Ycopen file for results 

fprintf(out, 'bi_dirty2 corr tests\n’); 

for testcase = 1:45 
% find most likely aircraft 
index=find(class_array {testcase} ==max(class_array {testcase} )); 
% fill matrix with coefficients 
coef_matrix_d2(rowD(testcase),:) = coef_matrix_d2(rowD(testcase),:) + fillRowWithCoefs(class_array {testcase} ); 
% fill binary matrix with likely identification 
bin_matrix_d2(rowD(testcase),column(index))=bin_matrix_d2(rowD(testcase),column(index))+1; 
fprintf(out, '%3g %s\n',testcase, lookup_table(index,:)); 

end 

clear R_dirty2; 

coef_matrix_d2 = coef_matrix_d2 / 5; %average matrix (5 cases per aircraft) 

fclose(out); 
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% Consolidate and print results to a file 


% results will be matrices where rows correspond to test case platforms (5 
% each) and columns correspond to library identifications by aircraft type. 


out = fopen(‘bi_matrix_bin.txt','wt'); 
fprintf(out, "bi_stock corr tests\n\n'); 


fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 
fprintf(out,'%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f 


fprintf(out, '\nbi_stock2 corr tests\n\n'); 


fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 
fprintf(out,'%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f 


fprintf(out, '\nbi_dirty corr tests\n\n'); 


fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 
fprintf(out, '%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f 


fprintf(out, '\nbi_dirty2 corr tests\n\n'); 





fprintf(out, 'A320 EA6B F4 F14 F15 Fl6 FA18 MiG Su27 UH60\n'); 
fprintf(out,'%1.0f %1.0f %1.0f %1.0f %1.0f %1.0f %1.0f 


felose(out); 


% Consolidate and print results to a file 
out = fopen(‘bi_matrix_coef.txt’,'wt'); 
fprintf(out, "bi_stock corr tests\n\n'); 


fprintf(out,'A320 EA6B F4 F114 FI5 F16 FAI18 MiG Su27 U 
fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5 


fprintf(out, '\nbi_stock2 corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F14 FI5 F116 FAI8 MiG 


fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %S5. 


fprintf(out, '\nbi_dirty corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F114 FI5 F116 FAI8 MiG 


fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5 


fprintf(out, '\nbi_dirty2 corr tests\n\n'); 
fprintf(out,'A320 EA6B F4 F114 FI5 F116 FAI8 MiG 














%1.0f %1.0f %1.0f \n', bin_matrix_s'); 


%LOF %1.0f %1.0f \n', bin_matrix_d’); 


%1.0f %1.0f %1.0f \n', bin_matrix_d2') 


H60\n'); 


Af\n', coef_matrix_s'); 


H60\n’); 
4f\n', coef_matrix_s2'); 


H60\n’); 


Af\n', coef_matrix_d'); 





H60\n’); 


fprintf(out, '%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n', coef_matrix_d2'); 


felose(out); 


% AR 2A is Ag 2S 2s 2g 2S 2 2g 26 2s 2k 2g oi 2k 2 2K ok oo checkLibraryCORRBi.m 3A os fg 24S AS os 2g 2g ag 2g 2s 2k 2c ois os 2k ois 2s 2k 2 ok 2 ok 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% RUN TIME: 10 Hours 


% checkLibrary is a function that provides an array where each element 


% indicates how similiar the signal is to the specific library signal. 


% The array element corresponds to a lookup table of all the library 
% signals. It must be done one part of the library at a time since loading 


% the entire library would require 8GB of RAM. 


% R-tests is now a cell array of the 2D Fourier transforms of the test 
% cases with num being the number of test cases contained in the the cell 


% array. Pre-transformed for speed. 


% Function flow loads the bispectrum array of 208 profiles for each 

% aircraft one by one and compares each of the library profiles with all 
% test cases. output_array becomes a cell array where each cell is an 
% array of correlation coefficients for that test case against all library 
% instances. The fft2 is not fftshifted because the result is not effected 


% by the change. 
% each aircraft_bic is over 700MB 
function output_array = checkLibraryCORRBi(R_tests,num); 


load ../zLibrary/a320_bic; %load a320 section of library 
for index=1:208 %for the a320 section of the output_array 
Rlib = fft2(abs(a320_bic {index})); 
% compare this library case to all test cases 
for test=1:num 
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%1.0f %1L.0f %1.0f \n', bin_matrix_s2'); 


> 


output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear a320_bic; 


load ../zLibrary/ea6b_bic; “load ea6b section of library 
for index=209:416 %for the ea6b section of the output_array 
Rlib = fft2(abs(ea6b_bic {index-208})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear ea6b_bic; 


load ../zLibrary/f4_bic; Yoload f4 section of library 
for index=417:624 %for the f4 section of the output_array 
Rlib = fft2(abs(f4_bic {index-416})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear f4 bic; 


load ../zLibrary/f14_bic; “load f14 section of library 
for index=625:832 %%for the f14 section of the output_array 
Rlib = fft2(abs(f14_bic {index-624})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear f14_bic; 


load ../zLibrary/f15_bic; “load f15 section of library 
for index=833:1040 %for the f15 section of the output_array 
Rlib = fft2(abs(f15_bic {index-832})); 
Y%compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear f15_bic; 


load ../zLibrary/f16_bic; “load f16 section of library 
for index=1041:1248 %for the f16 section of the output_array 
Rlib = fft2(abs(f16_bic {index-1040})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear f16_bic; 


load ../zLibrary/fal8_ bic; “load fal8 section of library 
for index=1249:1456 %for the fal8 section of the output_array 
Rlib = fft2(abs(fal8_bic {index-1248})); 
“% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear fal8_bic; 


load ../zLibrary/mig29_ bic; %load mig29 section of library 
for index=1457:1664 %for the mig29 section of the output_array 
Rlib = fft2(abs(mig29_bic {index-1456})); 


% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear mig29 bic; 


load ../zLibrary/su27_bic; Y%oload su27 section of library 
for index=1665:1872 %for the su27 section of the output_array 
Rlib = fft2(abs(su27_bic {index-1664})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear su27_bic; 


load ../zLibrary/uh60_bic; Yoload uh60 section of library 
for index=1873:2080 %for the uh60 section of the output_array 
Rlib = fft2(abs(uh60_bic {index-1872})); 
% compare this library case to all test cases 
for test=1:num 
output_array {test} (index) = howCorrBi(R_tests {test} ,Rlib); 
end 
end 
clear uh60_bic; 


% AS 2g os he 2s 2s he 24g 2 2g 2c os ok 2k os ok 2k 2 ok ok howCorrBi.m OAS hg 2S 2S is 2g 2A as 2g 2k 2s 2g 2k 2 2k 2k ok 2 2k 2k oi 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% howCorrBi determines how closely the test signal matches the library 
% signal. Inputs are 2D Fourier transforms of bispectra 


% input: test bicorr <511x511> and library bicorr <511x511> 
% output: corr coef between 0 and | 


function result = howCorrBi(Rtest,Rlib); 

% implements equation (15) from chapter V 

denominator = sqrt(sum(sum(Rtest.*conj(Rtest))))*sqrt(sum(sum(Rlib.*conj(Rlib)))); 
ccoef = ifft2(Rtest.*conj(Rlib))/denominator; 


result = abs(max(max(ccoef)))*511*511; 


% AR 2g ig fg 2A 2s 2g 2g 2 2g 2g 2s ok 2g os 2 2 oo oo fillRow WithCoefs.m OAS Sg 2S 2S fg 2g 2S fe 2g 2 os 2k 2g os 2k ois ok 2k 2 ok 2 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% this puts the max correlation coefficient for each 
% aircraft type into the appropriate place in the matrix row. 


function result = fillRowWithCoefs(arr); 
result = zeros(1,10); 


if length(arr) == 2080 “if stock cases 
result(1) = max(arr(1:208)); max of a320 checks 
result(2) = max(arr(209:416)); Yomax of ea6b checks 
result(3) = max(arr(417:624)); Yomax of £4 checks 
result(4) = max(arr(625:832)); Yomax of f14 checks 
result(5) = max(arr(833:1040)); Yomax of f15 checks 
result(6) = max(arr(1041:1248)); Yomax of f16 checks 
result(7) = max(arr(1249:1456)); Yomax of fal8 checks 
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result(8) = max(arr(1457:1664)); Y%omax of mig29 checks 

result(9) = max(arr(1665:1872)); Yomax of su27 checks 

result(10) = max(arr(1873:2080)); Yomax of uh60 checks 
end 


if length(arr) == 1872 “if dirty cases 
result(2) = max(arr(1:208)); Yomax of ea6b checks 
result(3) = max(arr(209:416)); Yomax of f4 checks 
result(4) = max(arr(417:624)); Yomax of f14 checks 
result(5) = max(arr(625:832)); Yomax of f15 checks 
result(6) = max(arr(833:1040)); Yomax of f16 checks 
result(7) = max(arr(1041:1248)); Yomax of fal8 checks 
result(8) = max(arr(1249:1456)); Yomax of mig29 checks 
result(9) = max(arr(1457:1664)); Yomax of su27 checks 
result(10) = max(arr(1665:1872)); Yomax of uh60 checks 

end 


% Ag 2A os he 2k 2 he fg 2s is 2g 2 os 2k oi os 2k oi ok ok ok column.m ORG sg 26 2S is 2g 2k os 2g 2k 2s ik 2s os ok 2k oie oi 2k oo ok 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% simple function to convert the index to a column for the confusion matrix 
% the index corresponds to a specific aircraft's area of the array. 
% result is the number of the library aircraft. 


function result = column(index); 


if index <= 208 
result = 1; %a320 
elseif index <= 416 
result = 2; %ea6b 
elseif index <= 624 
result = 3; %f4 
elseif index <= 832 
result = 4; %f14 
elseif index <= 1040 
result = 5; %f15 
elseif index <= 1248 
result = 6; %f16 
elseif index <= 1456 
result = 7; %fal8 
elseif index <= 1664 
result = 8; Y%omig29 
elseif index <= 1872 
result = 9; %su27 
else 
result = 10; %uh60 
end 


% AR 2g os Ag 2s 2 fe 2g 2s fe Ag 2s 2k 26 2 ok 2g oo 2 rowS.m OAS AS ig Ag 2S 2S hg 2S 2S ig 26 2s ik 2g 2 ok 2k 2g 2 2 ook 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% converts the testcase variable to a row on the column matrix for the 
% stock test cases 


function result = rowS(test); 


if test <= 5 
result = 1; %a320 
elseif test <= 10 
result = 2; %ea6b 
elseif test <= 15 
result = 3; %f4 
elseif test <= 20 
result = 4; %f14 
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elseif test <= 25 
result = 5; %f15 
elseif test <= 30 
result = 6; %f16 
elseif test <= 35 
result = 7; %fal8 
elseif test <= 40 
result = 8; Y%omig29 
elseif test <= 45 
result = 9; %su27 
else 
result = 10; %uh60 
end 


% AG Ag he 2g 2s fe 2g 2s fe 2g 2s ok 2g 2 ok 2k oi ok 2 ok rowD.m BAS 2g 2S 2s 2g 26 2s ig 2s 2s 2k 2c ois os 2k ois ok 2k 2 ok 2 


% LT Zachary Cole, USN 
% Naval Postgraduate School - Physics 


% converts the testcase variable to a row on the column matrix for the 
% dirty test cases 


function result = rowD(test); 


if test <= 5 
result = 2; %ea6b 
elseif test <= 10 
result = 3; %f4 
elseif test <= 15 
result = 4; %f14 
elseif test <= 20 
result = 5; %f15 
elseif test <= 25 
result = 6; %f16 
elseif test <= 30 
result = 7; %fal8 
elseif test <= 35 
result = 8; Y%omig29 
elseif test <= 40 
result = 9; %su27 
else 
result = 10; %uh60 
end 


E. INDIVIDUAL SIMULATION RESULTS 


High Bandwidth Results: 
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